include("geometry.jl")
include("misc.jl")
include("rhs.jl")

function main(p)
    # order of interpolate polynomials used for approximation
    P = p

    # set the simulation type, only support InsentroicVortex test now
    SimType = "InsentroicVortex"

    if SimType == "InsentroicVortex"
        MeshFile = "vortex.neu"
    else
        print("Simulation case using unknown")
        exit()
    end

    # set element-wise operators
    Op = Operator(P)

    # set geomatric informations
    Geo = Geometry(MeshFile, Op)

    # set compute parameters
    γ = 1.4; t_0 = 0; t_f = 10.0; tstep = 1

    # initilize the solutions
    u = Initialize(Op.Np, Geo.K, Geo.SPs, t_0, SimType)

    # compute initial time-step
    dt = EulerDt(P, u, Geo.vmapM, γ, Geo.Fscale)

    # define the  RK coefficents
    rk4a = [0.0,
            -567301805773.0/1357537059087.0,
            -2404267990393.0/2016746695238.0,
            -3550918686646.0/2091501179385.0,
            -1275806237668.0/842570457699.0]

    rk4b = [1432997174477.0/9575080441755.0,
            5161836677717.0/13612068292357.0,
            1720146321549.0/2090206949498.0,
            3134564353537.0/4481467310338.0,
            2277821191437.0/14882151754819.0]

    rk4c = [0.0,
            1432997174477.0/9575080441755.0,
            2526269341429.0/6820363962896.0,
            2006345519317.0/3224310063776.0,
            2802321613138.0/2924317926251.0]

    # storage for low storage RK time stepping
    rhsu = zeros(Float64, Op.Np, Geo.K, 4)
    resu = zeros(Float64, Op.Np, Geo.K, 4)

    # temporal integration
    t = t_0
    while t < t_f
        # Check if we need to adjust for final time step
        if t + dt > t_f
            dt = t_f - t
        end

        # SSK loop
        for i = 1:5
            # compute right hand side of compressible Euler equations
            rhsu  = RHS(u, Op, Geo, γ, t, SimType)

            # initiate and increment Runge-Kutta residuals
            resu = rk4a[i]*resu + dt*rhsu

            # update fields
            u = u + rk4b[i]*resu
        end

        # update time and compute new timestep
        t = t + dt; tstep = tstep + 1
        dt = EulerDt(P, u, Geo.vmapM, γ, Geo.Fscale)
    end
    ue = Initialize(Op.Np, Geo.K, Geo.SPs, t_f, SimType)

    # evaluate L2 norms
    ϵ = (u - ue).^2

    ρ_e = sqrt(sum(ϵ[:,:,1]))/(Op.Np*Geo.K)
    ρu_e = sqrt(sum(ϵ[:,:,2]))/(Op.Np*Geo.K)
    ρv_e = sqrt(sum(ϵ[:,:,3]))/(Op.Np*Geo.K)
    E_e = sqrt(sum(ϵ[:,:,4]))/(Op.Np*Geo.K)

    error = (ρ_e, ρu_e, ρv_e, E_e)
    u, ue, error
end